Format of the data

To analyze an experiment for batch effects, three tables need to be loaded into environment:

The package typically requires three datasets: 1) measurement table (data matrix), 2) sample annotation and 3) feature annotation table (optional). If you are familiar with Biobase package, these correspond to 1) assayData, 2) joined phenoData and protocolData, 3)featureData.

• Measurement table (data matrix), either wide or long format data frame. In the wide (matrix) format, referred in this vignette as data_matrix, rows represent features (for proteomics, peptides/fragments) and columns represent samples. In the long format, referred in this vignette as df_long, each row is a measurement of a specific feature (peptide, fragment) in the specific sample, requires at least three columns: feature ID column, measurement (intensity) column and sample ID column.

In this vignette, the essential columns have the following names:

feature_id_col = 'peptide_group_label'
measure_col = 'Intensity'
sample_id_col = 'FullRunName'
essential_columns = c(feature_id_col, measure_col, sample_id_col)

The names of the columns can be technology-specific. These column names are OpenSWATH tsv output format.

In the package, we provide the function to convert from long to matrix format (see section XXX “Utility functions”).

Note that sample IDs (column names in data_matrix) should match the values in sample ID column in sample_annotation and feature ID column values should match feature annotation table (here - peptide_annotation). For OpenSWATH tsv file, which is a long format data frame, peptide_annotation can be generated in the beginning of the analysis.

• Sample annotation is a data frame, where one row corresponds to one sample (run/file), and the columns contain information on biological and technical factors. Minimally, sample annotation has to contain sample ID column, at least one technical and one biological factor column, and biological ID column (unique ID for the biological replicate, which is repeated for each technical replicate). In our example data, these columns are:

1. sample_id_col = 'FullRunName'
2. technical covariates:

1. SacrificeDate - date when tissues were extracted
2. ProteinPrepDate - date when samples were prepared
3. RunDate (and RunTime, if available) - will be used to determine run order
4. MS_batch - number of MS batches (in this case, sets of runs between machine cleaning)

3. biological covatiates:

1. Strain
2. Diet
3. Sex

4. biospecimen_id_col = "EarTag" ’ Thus, technical and biological factors are:

technical_covariates = c('MS_batch', 'digestion_batch', 'RunDate', 'RunTime')
biological_covariates = c('Strain', 'Diet', 'Sex', 'Age_Days')
biospecimen_id_col = "EarTag"

In the example dataset, you will also find order column, which is, however, best inferred from the run date/time with the corresponding function (see Utility functions below).

• Feature (peptide) annotation: is a dataframe, where one row corresponds to one feature (in MS proteomics - peptide or fragment), and the columns are names of proteins and corresponding genes. Thus, the minimum is feature ID (peptide_group_id) and name of corresponding protein (in this vignette, we use Gene name).

Example dataset

The proBatch package can be applied to any dataset, for which intensity matrix and sample annotation are available, althogh the package was primarily designed with proteomic data in mind, and thoroughly tested on SWATH data. Thus, as and example data we include a reduced SWATH-MS measurement file, generated from BXD mouse aging study. In this study, liver proteome of mouse from BXD reference population have been profiled to identify proteome changes associated with age. The animals of each strain were subjected to Chow and High-Fat Diet and sacrificed at different age (age factor is excluded from the example data as age-related differences are in the focus of unpublished manuscript).

This dataset has a few features, that make it a good illustrative example: 1. This is a large dataset of 371 samples, that was affected by multiple technical factors, described above in sample_annotation subsection. Specifically, 7 MS batches drive the similarity of the samples. 2. The technical factors bias the data in at least two ways: discrete shifts affect different peptides in a batch-specific way, and MS drift, that introduces continous bias associated with sample running order. We will illustrate, how such biases can be corrected. 3. Replicate structure: two animals were injected every 10-15 samples. Additionally, several samples were repeated back-to-back in the end and in the beginning of two consecutive batches. This replication scheme allows to evaluate coefficient of variation and is highly beneficial for assessment of sample correlation.

The example SWATH data and annotation files can be loaded from the package with the function data().

library(proBatch)
data("example_proteome", "example_sample_annotation", "example_peptide_annotation", 
     package = "proBatch")

Prepare sample and peptide annotations

Before you skip this section, read whether you need the enhanced functionality that proBatch provides for preparing the annotation of samples and peptides.

We provide a few utility functions to facilitate the preparation of sample and peptide annotation.

generated_sample_annotation <- date_to_sample_order(example_sample_annotation,
                                          time_column = c('RunDate','RunTime'),
                                          new_time_column = 'generated_DateTime',
                                          dateTimeFormat = c("%b_%d", "%H:%M:%S"),
                                          new_order_col = 'generated_order',
                                          instrument_col = NULL)
library(knitr)
kable(generated_sample_annotation[1:5,] %>%
  select(c("RunDate", "RunTime", "order", "generated_DateTime", "generated_order")))

generated_peptide_annotation <- create_peptide_annotation(example_proteome, 
                                        feature_id_col = 'peptide_group_label',
                                        annotation_col = c('ProteinName', 'Gene'))

example_proteome = example_proteome %>% select(one_of(essential_columns))
gc()
#>           used  (Mb) gc trigger  (Mb) max used  (Mb)
#> Ncells 3168075 169.2    4703850 251.3  4703850 251.3
#> Vcells 3969774  30.3   14349595 109.5 16602139 126.7

Other utility functions

example_matrix <- long_to_matrix(example_proteome)

log_transformed_matrix <- log_transform(example_matrix)

color_scheme <- sample_annotation_to_colors (example_sample_annotation, 
       factor_columns = c('MS_batch','EarTag', "Strain", "Diet", "digestion_batch", "Sex"),
       not_factor_columns = 'DateTime',
       numeric_columns = c('Age_Days', 'order'))
color_list = color_scheme$list_of_colors

Plot sample mean

The plot_sample_mean illustrates global average vs. sample order. This can be helpful to visualize the global quantitative pattern and identify discrepancies within or between batches.

batch_col = 'MS_batch'
plot_sample_mean(log_transformed_matrix, example_sample_annotation, order_col = 'order', 
                 batch_col = batch_col, color_by_batch = T, ylimits = c(12, 16),
                 color_scheme = color_list[[batch_col]])

We can clearly see down-sloping trends in the BXD aging dataset. In fact, during the data acquisition, the mass-spectrometer had to be interrupted several times for tuning and/or column exchange as the signal was decreasing.

Plot boxplots

Alternatively, plot_boxplots() captures the global distribution vs. the sample running order.

log_transformed_long <- matrix_to_long(log_transformed_matrix)
batch_col = 'MS_batch'
plot_boxplot(log_transformed_long, example_sample_annotation, 
             batch_col = batch_col, color_scheme = color_list[[batch_col]])

In many cases, global quantitative properties such as sample medians or standard deviations won’t match. The initial assessment via mean plots or boxplots can capture such information and hint at which normalization method is better suitable. If the distributions are comparable, methods as simple as global median centering can fix the signal shift, while quantile normalization can help in case of divergent distributions.

Median normalization

Median normalization is a conservative approach that shifts the intensity of the sample to the global median of the experiment. If distributions of samples are dramatically different and this cannot be explained by non-technical factors, such as heterogeneity of samples, other approaches, such as quantile normalization need to be used.

median_normalized_matrix = normalize_data(log_transformed_matrix, 
                                          normalizeFunc = "medianCentering")

Same result will be achieved with:

median_normalized_matrix = normalize_data(example_matrix, 
                                          normalizeFunc = "medianCentering", log_base = 2)

Quantile normalization

Quantile normalization sets different distributions of individual samples to the same quantiles, which forces the distribution of the raw signal intensities to be the same in all samples. This method is computationally effective and has simple assumption that the majority of features (genes, proteins) is constant among the samples, thus also the distribution in principle are identical.

quantile_normalized_matrix = normalize_data(log_transformed_matrix, 
                                            normalizeFunc = "quantile")

After quantile or median normalization, you can easily check if the global pattern improved by generating mean or boxplots and comparing them side by side. Here are the mean plots before and after normalization of the log transformed dataset.

plot_sample_mean(quantile_normalized_matrix, example_sample_annotation, 
                 color_by_batch = T, ylimits = c(12, 16), 
                 color_scheme = color_list[[batch_col]])

Hierarchical clustering

Hierarchical clustering is an algorithm that groups similar samples into a tree-like structure called dendrogram. Similar samples cluster together and the driving force of this similarity can be visualized by coloring the leaves of the dendrogram by technical and biological variables.

Our package provides plot_sample_clustering() and plot_heatmap() to plot the dendrogram by itself or with a heatmap. You can easily color annotations on the leaves of the dendrograms or heatmaps to identify what is the driving force of clustering.

Once your color annotation is ready, for the specific covariates of interest, you can subset the color dataset and feed it into the clustering functions.

library(dplyr)
color_annotation <- color_scheme$color_df
selected_annotations <- c("MS_batch",  "digestion_batch", "Diet")
#select only a subset of samples for plotting
color_annotation <- color_annotation %>% select(one_of(selected_annotations))

#Plot clustering between samples 
plot_hierarchical_clustering(quantile_normalized_matrix, color_annotation,  
                             distance = "euclidean", agglomeration = 'complete',
                             label_samples = F)

Similarly, you can plot a heatmap by supplementing the color list. You decide whether to show annotations in the column, row or both by specifying required covariates in sample_annotation_col, sample_annoation_row, or both.

plot_heatmap(quantile_normalized_matrix, example_sample_annotation, 
             sample_annotation_col = selected_annotations, 
             cluster_cols = T, 
             annotation_color_list = color_scheme$list_of_colors,
             show_rownames = F, show_colnames = F)

From the clustering analysis, we can clearly see that the driving force behind sample clustering is the MS batch.

Principal component analysis (PCA)

PCA is a technique that identifies the leading directions of variation, known as principal components. The projection of data on two principal components allows to visualize sample proximity. This technique is particularly convenient to assess replicate similarity.

You can identify the covariate leading the direction of variations by coloring potential candidates.

plot_PCA(quantile_normalized_matrix, example_sample_annotation, color_by = 'MS_batch', 
              plot_title = "MS batch", colors_for_factor = color_list[[batch_col]])
plot_PCA(quantile_normalized_matrix, example_sample_annotation, color_by = "digestion_batch", 
         plot_title = "Digestion batch", colors_for_factor = color_list[["digestion_batch"]])
plot_PCA(quantile_normalized_matrix, example_sample_annotation, color_by = "Diet",  
         plot_title = "Diet", colors_for_factor = color_list[["Diet"]])

By plotting the first two principal components and applying different color overlaps, we see once again that clusters overlap nicely with MS batches.

Principal variance component analysis (PVCA)

The main advantage of this approach is the quantification of the variance, associated with both technical and biological covariates. Briefly, principal variance component analysis uses a linear model to match each principal component to the sources of variation and weighs the variance of each covariate by the eigenvalue of the PC [19]. Thus, the resulting value reflects the variance explained by that covariate.

NB: PVCA calculation is a computationally demanding procedure. For the data matrix of several hundred samples and several thousands of peptides it can easily take several hours. So if it is generally a good idea to run this analysis as a stand-alone script on a powerful machine.

pvca <- plot_PVCA(quantile_normalized_matrix, example_sample_annotation, 
                  technical_covariates = c('MS_batch', 'digestion_batch'),
                  biological_covariates = c(biological_covariates, biospecimen_id_col))

The biggest proportion of variance in peptide measurement was derived from mass spectrometry batches. In typical experiment, the overall magnitude of variances coming from biological factors should be high while technical variance should be kept at minimum¹.

Peptide-level diagnostics and spike-ins

Featre-level diagnostics is very informative for batch effect correction. To assess the bias in the data, one can choose a feature (peptide, protein, gene), quantitative behaviour of which is known. Through our package, plot_peptides_of_one_protein() allow plotting peptides of interest e.g. a biologically well understood protein. If spike-in proteins or peptides have been added to the mixture, one can use the plot_spike_ins() function. In most DIA datasets iRT peptides [22] are added in controlled quantities and can be visualised with plot_iRT() function.

In mass-spectrometry, also the trends associated with this order can be assessed for a few representative peptides, thus order column is also important for this diagnostics.

quantile_normalized_long <- matrix_to_long(quantile_normalized_matrix, example_sample_annotation)
plot_spike_in(quantile_normalized_long, example_sample_annotation, 
              peptide_annotation = generated_peptide_annotation,
              protein_col = 'Gene', spike_ins = "BOVINE_A1ag", 
              plot_title = 'Spike-in BOVINE protein peptides',
              color_by_batch = T, color_scheme = color_list[[batch_col]])

It is clear that while the pre-determined quantities of spike-ins or peptides of known biology have their expected intensities, the trend is dominated by mass spectrometry signal drift. After confirming either continuous or discrete batch effects exist in a dataset, by one or more of these methods, proceed by selecting a batch correction method.

Continuous drift correction

Continuous drifts are specific for mass-spectrometry, thus, the user-friendly methods for its correction have not been implemented before. In this package, we suggest a novel procedure to correct for MS signal drift. We developed a procedure based on nonlinear LOESS fitting. For each peptide and each batch, a non-linear trend is fitted to the normalized data and this trend is subtracted to correct for within-batch variation. Note that the resulting data are not batch-free as within-batch means and variances are batch-dependent. However, now the batches are discrete and thus can be corrected discrete methods such as median-centering or ComBat.

loess_fit <- adjust_batch_trend(quantile_normalized_matrix, example_sample_annotation)
loess_fit_matrix <- loess_fit$data_matrix

One important parameter in LOESS fitting is span, which determinesthe degree of smoothing. The LOESS span ranges from 0 to 1: the greater the value of span, the smoother is the fitted curve. Since we want the curve to reflect signal drift, we want to avoid overfitting but be sensitive to fit the trend accurately. Currently, we suggest to evaluate several peptides to determine the best smoothing degree for a given dataset.

loess_fit_30 <- adjust_batch_trend(quantile_normalized_matrix, example_sample_annotation, span = 0.3)

quantile_normalized_long <- matrix_to_long(quantile_normalized_matrix)
plot_with_fitting_curve(pep_name = "10231_QDVDVWLWQQEGSSK_2", 
            df_long = quantile_normalized_long, example_sample_annotation, 
            color_by_batch = T, color_scheme = color_list[[batch_col]],
            fit_df = loess_fit_30$fit_df, plot_title = "Span = 30%")

loess_fit_70 <- adjust_batch_trend(quantile_normalized_matrix, example_sample_annotation, span = 0.7)
plot_with_fitting_curve(pep_name = "10231_QDVDVWLWQQEGSSK_2", 
            df_long = quantile_normalized_long, example_sample_annotation, 
            color_by_batch = T, color_scheme = color_list[[batch_col]],
            fit_df = loess_fit_70$fit_df, plot_title = "Span = 70%")

Curve fitting is largely dependent on number of consecutive measurements. In proteomics, missing values are not uncommon. If too many points are missing, a curve cannot be fit accurately. This is especially common for small batches. In this case, we suggest to not fit the curve to the specific peptide within the specific batch, and proceed directly to discrete correction methods. To identify such peptdies, absolute and relative threshold (abs.threshold and pct.threshold) on number of missing values can be substitutes as paramters for normalize_custom_fit() is the absolute and relative threshold on the number of missing values for each peptide.

Discrete batch correction: combat or peptide-level median centering

Once the data are normalized and corrected for continuous drift, only discrete batch effect is left to be corrected.

Currently, both batch correction procedures accept long format of the data:

loess_fit_df <- matrix_to_long(loess_fit_matrix)

peptide_median_df <- center_peptide_batch_medians(loess_fit_df, example_sample_annotation)
plot_single_feature(pep_name = "46213_NVGVSFYADKPEVTQEQK_2", df_long = peptide_median_df, 
            example_sample_annotation, color_by_col = NULL, measure_col = 'Intensity_normalized',
            plot_title = "Feature-level Median Centered")

comBat_matrix <- correct_with_ComBat(loess_fit_matrix, example_sample_annotation)
#> Standardizing Data across genes

combat_df <- matrix_to_long(comBat_matrix)
plot_single_feature (pep_name = "46213_NVGVSFYADKPEVTQEQK_2", loess_fit_df, 
          example_sample_annotation, plot_title = "Loess Fitted", color_by_col = NULL)
plot_single_feature (pep_name = "46213_NVGVSFYADKPEVTQEQK_2", combat_df, 
          example_sample_annotation, plot_title = "ComBat corrected", color_by_col = NULL)

Correct batch effects: universal function

We provide a convenient all-in-one function for batch correction. The function correct_batch_effects corrects MS signal drift and discrete shift in a single function call. Simply specify which discrete correction method is preferred at discreteFunc either “ComBat” or “MedianCentering” and supplement other arguments such as span, abs.threshold or pct.threshold as in normalize_custom_fit().

batch_corrected_matrix <- correct_batch_effects(data_matrix = quantile_normalized_matrix, 
                                  example_sample_annotation, discreteFunc = 'ComBat',
                                  abs.threshold = 5, pct.threshold = 0.20)
#> Standardizing Data across genes

Correlation distribution of samples

For example, in the mice aging experiment, biological replicates, ET1506 and ET1524, were injected every 30-40 MS runs. The correlation between these biological replicates should improve after normalization and batch correction.

The plot_sample_corr_distribution() plots correlation distribution between biological replicates and non-replicates in the same or differnet batches by plot_param = 'batch_replicate'. Alternatively, you can compute the correlation between different batches by plot_param = 'batches'.

It should be noted, however, that the comparison of sample correlation should be approached with the comparison of distributions, rather than by evaluating the individual examples of within-replicate vs within-batch corrections. The sample correlation is often used to prove the quality of the measurement, as it is typically very high (examples of the replicate correlation above .95 are common for mass spectrometry). However, unless these examples are shown in the context of the whole distribution structure, they can lead to erroneous conclusion.

sample_cor_norm <- plot_sample_corr_distribution(quantile_normalized_matrix,
                                                 example_sample_annotation, 
                                                 batch_col = 'MS_batch', 
                                                 biospecimen_id_col = "EarTag", 
                                                 plot_title = 'Quantile normalized',
                                                 plot_param = 'batch_replicate')
sample_cor_norm + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + ylim(0.7,1)

sample_cor_batchCor <- plot_sample_corr_distribution(batch_corrected_matrix,
                                                     example_sample_annotation, 
                                                     batch_col = 'MS_batch', 
                                                     plot_title = 'Batch corrected',
                                                     plot_param = 'batch_replicate')
sample_cor_batchCor + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + ylim(0.7, 1)

Correlation distribution of peptides within and between proteins

Especially for proteins of biological significance, peptides of the same protein should behave similarly in expression level. Therefore, we can compare peptide correlation before and after batch correction if the correlation between the same proteins increases while that of different proteins stays the same.

peptide_cor_norm <- plot_peptide_corr_distribution(quantile_normalized_matrix, 
          generated_peptide_annotation, protein_col = 'Gene', plot_title = 'Quantile normalized')
peptide_cor_norm + geom_hline(yintercept=0, linetype="dashed", color = "grey")

peptide_cor_batchCor <- plot_peptide_corr_distribution(batch_corrected_matrix, 
          generated_peptide_annotation, protein_col = 'Gene', plot_title = 'Batch corrected')
peptide_cor_batchCor + geom_hline(yintercept=0, linetype="dashed", color = "grey")

Heatmap of selected replicate samples

In this study, 10 samples were run in the same order before and after the tuning of the mass-spectrometer, which marks the boundary between batches 2 and 3. The correlation between these replicates can be illustrated by corrplot (flavor = 'corrplot') or by pheatmap (flavor = 'pheatmap')

library(RColorBrewer)
earTags <- c("ET1524", "ET2078", "ET1322", "ET1566", "ET1354", "ET1420", "ET2154",
             "ET1515", "ET1506", "ET2577", "ET1681", "ET1585", "ET1518", "ET1906")

# sample ID of biological replicates 
bioreplicates <- example_sample_annotation %>%
  filter(MS_batch == 'Batch_2' | MS_batch == "Batch_3") %>%
  filter(EarTag %in% earTags) %>%
  pull(FullRunName)

# Prepare color annotation 
color_covariates <- example_sample_annotation %>%
  filter(MS_batch == 'Batch_2' | MS_batch == "Batch_3") %>%
  filter(FullRunName %in% bioreplicates) %>%
  remove_rownames %>% 
  column_to_rownames(var="FullRunName") %>%
  select("MS_batch", "EarTag") # Annotate MS_batch and EarTag on pheatmap 

# Plot pheatmap 
breaksList <- seq(0.7, 1, by = 0.01) # color scale of pheatmap 
plot_sample_heatmap(quantile_normalized_matrix, annotation_col = color_covariates, 
    samples_to_plot = bioreplicates, flavor = 'pheatmap', 
    color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(length(breaksList)),
    breaks = breaksList, plot_title = 'Quantile Normalized', cluster_rows=F, cluster_cols=F,
    annotation_names_col = TRUE, annotation_legend = TRUE, show_colnames = F)

plot_sample_heatmap(batch_corrected_matrix, annotation_col = color_covariates,
    samples_to_plot = bioreplicates, flavor = 'pheatmap', 
    color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(length(breaksList)),
    breaks = breaksList, plot_title = 'Batch Corrected', cluster_rows=F, cluster_cols=F,
    annotation_names_col = TRUE, annotation_legend = TRUE, show_colnames = F)

It is clear that before the correction, samples from one batch cluster better than between the replicates. However, after the correction, the correlation between replicates becomes higher than the correlation between non-related samples regardless of the batch.